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1. Introduction 



The top quark is the heaviest elementary particle known to date. Its large mass of 173 GeV 
makes it a crucial tool for the study of the electroweak symmetry breaking mechanism and 
for the search of signals of physics beyond the Standard Model. The study of the properties 
of the top quark at the Tevatron during the last 15 years provided precise measurements [1] 
of the top-quark mass (with a relative error of 0.75%) and of the total top-quark pair 
production cross section (with a relative error of about 10%). Owing to the by- now large 
Tevatron dataset, measurements of differential distributions in top pair production are 
becoming possible [2,3]. A new set of precise experimental measurements is expected 
to be obtained by the LHC experiments soon. At the LHC, running at 7 TeV, with an 
integrated luminosity of only ~ 200 pb -1 , it should be possible to record about 30000 top 
quark pair events before selection [4]. With the planned increases in center of mass energy 
and luminosity, the LHC will turn into a veritable top-quark factory, producing millions of 
top-quark events per year [5] . 

In order to take full advantage of the improved experimental measurements, it will 
be necessary to provide theoretical predictions for the measured observables which are as 
accurate as possible. Most of the observables related to the top-quark pair production have 
been calculated up to NLO [6]. In several cases, the next-to-leading (NLO) corrections have 
been supplemented by the resummation of large logarithmic corrections at leading (LL, 
[7]), next-to-leading (NLL, [8]) and next-to-next-to-leading logarithmic (NNLL, [9-12]) 
accuracy. However, to match the precision of the forthcoming experimental data, full 
next-to-next-to-leading order (NNLO) calculations are required for at least some of the 
observables, such as the top-quark pair production total cross section [13]. 

To obtain NNLO predictions for the top quark pair production in QCD, the following 
matrix elements need to be computed: i) two- loop matrix elements in both the quark- 
antiquark annihilation channel and the gluon-fusion channel; ii) one-loop matrix elements 
with one extra parton in the final state; and Hi) tree- level matrix elements with two extra 
partons in the final state. The diagrams belonging to the sets ii) and Hi) form part of 
the NLO corrections to top-pair-plus-jet production and have been calculated by several 
groups [14-16]. Aside from the calculation of the required Feynman diagrams, a full NNLO 
calculation of the top-quark pair production cross section is particularly challenging because 
it entails the development of a NNLO subtraction method, which could be used in processes 
with massive partons [17-21]. 

As far as matrix elements are concerned, the last missing ingredient to the NNLO 
corrections to top quark pair production are the two-loop matrix elements for qq — > ti 
and gg — > ti. Both types of matrix elements were calculated in the s S> m 2 limit (where 
s is the partonic center of mass energy and m is the mass of the top quark) [22, 23]. A 
fully numerical calculation of the two-loop corrections in the quark-antiquark annihilation 
channel is also available [24]. 

This paper is the third of a series of works aiming towards the analytic calculation 
of the two-loop matrix elements for the top pair production process. For what concerns 
the quark-antiquark annihilation channel, the two-loop diagrams involving a closed light or 
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heavy-quark loop were evaluated in [25], while the two-loop diagrams contributing to the 
leading color coefficient were evaluated in [26]. In both cases, the results obtained retain the 
full dependence on the top-quark mass and on the kinematic invariants; they agree with 
the numerical results of [24]. Having analytical results available has several advantages 
over a purely numerical representation. Besides their considerably shorter evaluation time, 
the analytical results also allow for an expansion in different kinematical limits (threshold, 
high energy). 

In the present paper, an analytical expression for the two-loop diagrams contribut- 
ing to the leading color coefficient in the gluon-fusion channel is derived. We carry out 
the calculation by employing the technique based on the Laporta algorithm [27] and the 
differential equation method [28], already used in [25,26]. The calculation of the leading 
color coefficient in the gluon fusion does not require the calculation of any new master 
integrals beyond the ones obtained in the two previous works, such that we do not discuss 
the calculational method in full detail. The interested reader can find in [25,26] a detailed 
description of the techniques employed. 

The paper is organized as follows: in Section ||, we introduce our notation and conven- 
tions; in Section [| we summarize the most relevant features of our calculational method. 
Section ||| describes the UV renormalization of the bare amplitude. The resulting two-loop 
amplitude contributions are described in Section || where we also provide numerical val- 
ues in some benchmark points, and discuss the expansion in the threshold limit. Finally, 
Section || contains our conclusions. 

2. Notation and Conventions 

We consider the scattering process 

g(pi) + g{P2) — ► t(p 3 ) + t(p 4 ) , (2.1) 

in Euclidean kinematics, where pf = for 2 = 1,2 and p 2 = —m 2 for j = 3,4. The 
Mandelstam variables are defined as follows 

s = ~ (jPi +P2) 2 , t = -(p 1 -p 3 ) 2 , u = -(pi-p i ) 2 . (2.2) 

Conservation of momentum implies that s + t + u = 2m 2 . 

The squared matrix element (summed over spin and color), calculated in d = 4 — 2e 
dimensions, can be expanded in powers of the strong coupling constant as as follows: 

\M\ 2 (s,t,m,e) = 16tt 2 q| 



Ai + 



7T 



A 2 + O (af) 



(2.3) 
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The tree-level amplitude involves the three diagrams shown in Fig. |T| and their contribution 
to Eq. ( p. 3| ) is given by 



Ao = (N 2 - 1 



1 1 



Nr. 



-N 



N c t 2 u 2 ° A u \ (*i + u l) 2 
4 [t\ + u\) (if + uiti + u?) 



Fi(ti,ni,m 2 ) 

l 2(ti+m) 2 



1 4 (if + uiti + u\) 



N r 



N c t\ui 
2{t\ + ul)- N 



iini(ti + ni) 2 J " L^c ii^i 

where iV c is the number of colors, t\ = t — m 2 , u\ = u — m 2 , and 

Fx(ti,ux,m) = -2txux(tx +uj) + 8m 2 t 1 u 1 (t l + ux) + 8m 4 (t 1 + m) 2 



(2.4) 



(2.5) 



The O(as) term Ax in Eq. fl2.3|) arises from the interference of one-loop diagrams with 
the tree-level amplitude [6]. The 0{aV) term Ai consists of two parts, the interference 
of two-loop diagrams with the Born amplitude and the interference of one-loop diagrams 
among themselves: 

A 2 



.(2X0) ,(lxl) 
•^2 ' "^2 



The term A^*^ was derived in [29-32]. A^*"' , originating from the two-loop diagrams 
can be decomposed according to color and flavor structures as follows: 

(N, 



(2x0) 



-(2x0) 



1)\n*A + N c B + + -^D + N 2 N t Ei + N 2 N h .L h + X, F, + V,, f), 



Nr 



iV 2 



N c N c 



(2.6) 



where N and iV^ are the number of light- and heavy-quark flavors, respectively. The coef- 
ficients A, B Iih in Eq. ( |2.6| ) are functions of s, t, and m, as well as of the dimensional 
regulator e. These quantities were calculated in [23] in the approximation s, \t\, \u\ ^> m 2 . 
For a fully differential description of top quark pair production at NNLO, the complete 
mass dependence of A^*^ is required; to date no such a result is available. Recently, a 
formula which allows to calculate the IR poles of a generic two-loop QCD amplitude was 
derived [33, 34] and employed to obtain analytic expressions for all of the IR poles for the 
top-quark pair production [35]. 

In this work, we provide an exact analytic expression for the leading colour coefficient 
A in Eq. (|2.6D, which receives contributions from planar Feynman diagrams only. 



3. Calculation 

The package QGRAF [36] generates 789 two-loop Feynman diagrams contributing to the 
process gg — > tt (considering one massless and one massive flavor). Evaluating their color 
structures, we find that 300 of them contribute to the leading color coefficient in Eq. fl2.G|). 
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Figure 2: Some of the two-loop planar box diagrams involved in the calculation. 

Since we use the covariant sum over the polarizations of the incoming gluons, we have 
to consider additional 116 (out of 218) diagrams for the process initiated by ghosts. We 
interfere the two-loop diagrams with the tree-level amplitude, and calculate the color and 
Dirac traces before carrying out the integration over the loop momenta. The resulting scalar 
loop integrals are reduced to a set of Master Integrals (Mis) employing the technique based 
upon the Laporta algorithm [27]. Then, the Mis can be evaluated analytically by means 
of the differential equation method [28]. 

Starting with the QGRAF output, the calculation of the interferences in terms of the 
Mis was carried out with a parallelized C++ package: Reduze 2 [37]. In particular, this 
code provides a fully distributed variant of the Laporta algorithm for the reduction 1 of the 
loop integrals. The package employs GiNaC [41] for the algebraic manipulations. 

The Mis needed for the calculation presented in this paper were already known in the 
literature [25,26,42-48]. In particular, all of the two-loop four-point Mis encountered in the 
calculation of the leading color coefficient in the gluon-fusion channel coincide with the ones 
needed for the corresponding calculation in the quark-antiquark annihilation channel [26], 
or can be obtained by the latter by replacing the Mandelstam variable t with u. All of the 
Mis were calculated in the non-physical region s < 0, where they are real. The result in 
the physical region is then recovered by analytic continuation. 

The transcendental functions appearing in the results are one- and two-dimensional 
harmonic polylogarithms (HPLs) [46, 49-51] of maximum weight four. All of the HPLs 
appearing in the analytic expression of the coefficient A can be evaluated numerically with 
arbitrary precision by employing the methods and codes described in [50]. A detailed list 
of the functional basis employed in this calculation can be found in Appendix A of [26]. 
Appendix B in the same paper describes the expansion of this class of HPLs in the thresh- 
old limit /3 — > 0, which requires some care. As expected, the sum of the bare two-loop 
corrections (as well as the UV renormalized corrections) is symmetric with respect to the 
exchange t o u. 



4. Renormalization 

The renormalized QCD amplitude is obtained from the bare one by expanding the following 

1 Other reduction codes released publicly are the Maple package A.I.R. [38], the Mathematica package 
FIRE [39] and the C++ package Reduze [40] . 
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expression : 

Aen = Y\ Z wf,u Aare ( a 5,bare ~^ Z a s a S , m bare ~^ Z m rnj , (4.1) 

n 

= Z WF ^ g Z WF j Abare ( Q S,bare ~^ Z a s ois ,^-bare ~* Z m rnj , (4.2) 

where Z WFjn (n = g,t) are the external particle wave-function renormalization factors, as 
is the renormalized coupling constant, Z as is the coupling constant renormalization factor, 
m is the renormalized heavy-quark mass, and Z m is the mass renormalization factor. In 
the rest of the section we suppress the subscript "S" in as- 
We introduce the following auxiliary quantities: 

TT TT 

By expanding the amplitude and the wave function renormalization factor in ao we find: 

Aen(«bare) = a oA + OqA + a^A 2 + 0{al) , 

Z WF , n = 1 + ao^w*F,n + Oo^WF.n + C(«o) > ( n = 9,t), 

Z m = l + aoSZW + a 2 ^ 2 ) + O(ag) . (4.4) 

The relation between ao and a is given by: 

a = a + a 2 5z2 ] + a 3 5Z^ + 0(a 4 ) . (4.5) 

By employing Eqs. (fO,4.5) in Eq. ( |4.2j ) we find 



Aen = aAo + a 2 A^ n + a 3 A ( $ n + C(a 4 ) , 

4'eli = A + (*zgj, t + + 25Z«) A + (^S lt ^ lfl + 2oZ«/Z« 

+25Z£IJZW + ^ + + 54 2 )) A - «#>.aJ~ CT) 

- (5Zi 2 ) + 2«5Z««5Z« + 6Z® tt 8Z£>) CT) + (^i 1 ^ A? CT) • (4-6) 

In the equations above, A> i = 0,1,2, represents the bare i-loop amplitude stripped of 
the factor a, .4,Q mass ^ represents the amplitude obtained by adding a mass insertion to 
the virtual top propagator in diagrams (b) and (c) of Fig. § A { ™ s ct) can be obtained 
by considering all possible mass insertions in the one- loop diagrams. Finally, Aq mass ^ is 
obtained by adding two mass insertions to the virtual propagator in diagrams (b) and (c) 
of Fig. |. 

In this work we employ a mixed renormalization scheme in which the wave functions 
and the heavy-quark mass are renormalized on shell, while the strong coupling constant is 
renormalized in the MS scheme. The explicit expressions of the one-loop renormalization 
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1/s 4 


1/e 3 


l/ £ 2 


1/s 


e° 


m 2 


m 2 


m 




5 


-1.25 


1 


10.749426 


18.693893 


-156.82372 


262.14826 


12.721807 


43 


-21 


1.7 


9.3642942 


-27.358589 


-41.372387 


305.80422 


-707.01281 


8.1 


-0.6 


2.1 


27.306074 


138.38466 


-381.11526 


-571.56293 


1859.6594 



Table 1: Numerical values of the coefficients appearing in the e expansion of A in the normalization 
that factors out the term [(47r) e e _7e ] . 



factors in Eq. ([4.6]) are: 



5Z, 



5Z 



(i) 

WF,t 
(1) 



WP,g — j 

^ = C(e 



-7£ 



T(l + e) V 2e 

0Z/ m — 0Z/ WF.t ' 



(4.7) 
(4.8) 
(4.9) 
(4.10) 



where C(e) = (4vr) £ r(l + e), /3 = (11/6)C A - (1/3) (JVj + JV h ) and where 7 is the Euler- 
Mascheroni constant 7 ~ 0.577216. 

Concerning the two-loop renormalization factors, we provide only the terms which 
contribute to the renormalization of the leading color coefficient. They are (see for in- 
stance [52]): 



5Z, 



(2) 



WF,i 



C(e? £j 



2 \ 2e 



C F 



( '1 ( i + ^~ + ^1 - ?C(3) + 6C(2) In 2 - f C(2) 



32e 2 64e 128 2 



(2) 
WF,g 



+C A 



0, 



11 £1 ™ » + « f(il)1 



32e 2 64e 128 4 



5Z& = C(e) 



-7£ 



1 

r(l + e)y 4e 

2 \ 2e 



6 



12 



(4.11) 

(4.12) 
(4.13) 



^ = c( £ ) 2 (4) ^ 



9 45 199 15 3 
C F [ -L- + — + — - — C(2) - -C(3) + 3C(2) In 2 
1 32e 2 64e 128 8 v ; 4 SV ; sv ; 



+C A 



(4.14) 



5. Results 



The main result of the present paper is an analytic expression for the coefficient A in 
Eq. ( |2.6| ) which retains the complete dependence on the variables t, u, on the renormaliza- 
tion scale /i, and on the top-quark mass m. The analytic result is too long to be explicitly 
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presented in written form. To make it available to the reader we include in the arXiv 
submission of this work a text file with the complete result. The coefficients in Eq. ( |2.C| ) 
still contain infrared poles. This makes our result dependent on the choice of a global, 
e-dependent normalization factor. We choose to present the coefficient A factoring out an 
overall term 

C 2 (e)= [(4^) £ r(l + £ )] 2 . (5.1) 

Another possible normalization, widely used in the literature, is the one in which a term 

C\e) = [(4vr) £ e -^] 2 (5.2) 

is retained as an overall factor (see for instance [35] and [24]). Our result, therefore, can 
be compared to the results in [24, 35] by multiplying it by a coefficient 

e 2 ^ T 2 (l + e) = 1 + C(2)e 2 - ^C(3)e 3 + ^C(2) V + O {e 5 ) . (5.3) 

While the formulae in the text and the supplied files are given in the normalization defined 
by Eq. ( |5.l|) , for ease of comparison with other results in the literature we present in Table [I] 
and in Figures a few benchmark points and plots in the normalization of Eq. (|5.2j). 

We provide a code which numerically evaluates the analytic expression of the leading 
color coefficient for arbitrary values of the mass scales involved in the calculation. The 
code is written in C++ and uses the package for the evaluation of multiple polylogarithms 
within GiNaC [50]. 

Our result was cross checked by comparing it with the partial results already available 
in the literature. By expanding the analytic expression of the coefficient A in the limit 
s, \t\, \u\ S> m 2 and neglecting terms suppressed by positive powers of m 2 /s for t/s = const, 
we obtained the result published in [23]. For the IR poles of the coefficient A we find 
numerical agreement with the results of [35]. 

In Fig. H the finite part of A is plotted as a surface depending on the variables n and 
(f), defined as 

*=A-i, * = J — , ifi-,/^W<i(i + ,/^V <M) 



4m 2 ' r s ' 2\ V 1 + r // 2 V V 1 + T ?, 

Finally, it is possible to expand the expression for A for values of the center of mass 
energy close to the production threshold. We define 



4m 2 1 — cos 9 ( fi 2 

— , L M = In — 

2 \m- 



where 9 is the scattering angle in the partonic center of mass frame. Keeping £ = const, we 
expand our results in powers of the heavy quark velocity (3 up to terms of order /3 5 included. 
The coefficients of this expansion contain transcendental constants which originate from 
one- and two-dimensional HPLs evaluated at x = 1. Since we did not find a satisfactory 
analytical representation for all of these constants, in the formulae below we provide the 
coefficients of the expansion in a numerical form. We find: 
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0.25 \ 





Figure 3: Finite part of the coefficient A as a function of the variables r\ = s/(4to 2 ) — 1 and 
(f) = — (t — m 2 )/s. We employed the normalization that factors out the term [(47r) e e _7£ ] and set 
fj, = m. 



with 

A { ~ A) = 8 + [24 - 32£(1 - £)] /3 2 + [32 - 128£(1 - £) - 256(£(1 - £)) 2 ] /3 4 

+0(/3 6 ) , (5.7) 

,4 ( - 3) = 13.1526 + 16L^ + [71.4579 + 48L M - (84.6105 + 64L^)£(1 - £)] P 2 
+ [103.277 + 64L /t - (487.775 + 256L /t )£(l - £) 

-(1850.22 + 512L M )(C(1 - £)) 2 )] /? 4 + ©(/3 6 ) , (5.8) 

A^ 2 ) = -155.954 - 32.3614L /t + 16L 2 

- [388.711 + 33.0843L M - 48L 2 - (1148.69 + 65.4457^ - 64L 2 )£(1 - £)] /3 2 

- [461.697 + 28.1123L M - 64L 2 - (2976.85 - 36.884L M - 256Lj)£(l - £) 
+(2814.61 + 1823.1L M + 512L 2 )(£(1 - £)) 2 ] /3 4 + 0(/3 6 ) , (5.9) 

A^-V = 224.144 - 191.469L. - 61.6948L 2 + 10.6667L?, 

t* 

+ [117.543 - 592.101L M - 121.084L 2 + 32L 3 

-(774.295 - 1463.61^ - 182.779L 2 + 42.6667L 3 )£(1 - £)] /3 2 

- [92.4819 + 744.745L M + 145.446L 2 - 42.6667L 3 
+(553.64 - 4574.2L^ - 432.449L 2 + 170.667L 3 )£(1 - f ) 
-(17463.5+2328.24L M -884.435iJ-341.333Lj)(^(l - £)) 2 ] /3 4 

+G(/3 6 ), (5.10) 
A (0) = 403.869 + 787.434L M - 77.4706L 2 - 50.9076L 3 + 5.33333L 4 

+ [1627.76 + 1320.97L M - 338.107L 2 - 110.056L 3 + 16L 4 

-(8363.19 + 4515.06L^ - 831.618L 2 - 160.964LJ + 21.3333L 4 )£(1 - £)] (3 2 

+ [1826.76 + 1164.29L M - 440.308L 2 - 136.075L 3 + 21.3333L 4 

-(22884.9 + 10052.4L^ - 3024.01L 2 - 444.744L 3 + 85.3333^)^(1 - £) 
+(50837.4+35185.&L, 4 +4586.08Lj-276.734Lj-170.667Lj)(^(l - £)) 2 ] /3 4 

+0(/3 6 ). (5.11) 
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Figure 4: Left: finite part of the coefficient A as a function of (3 = y/1 — 4m?/ s for £ = 1/2. 
Right: finite part of the coefficient A as a function of the variable r\ = s/(4m 2 ) — 1 for 0=1/2. In 
both cases we used the normalization that factors out the term [(4-7r) e e" 7£ ] and set fi = m. 



We observe that the dependence on j3 and on £ in the formulas above is polynomial. 
Moreover, [3 enters only via powers of /3 2 and £ only via powers of £(1 — £). The latter 
dependence explicitely reflects the symmetry of A under exchange of forward and backward 
directions, £ — > 1 — £. All of the logarithmic terms hi/3, ln£, ln(l — £), ln(l — 2^), . . ., 
which are indeed present in the expansion of individual HPLs, cancel out in the final 
expressions. Thus, the coefficient A is finite at threshold. The expansion presented here 
could be used in the future for the calculation of logarithmically enhanced leading colour 
terms near the tt production threshold. 

In Fig. ||, we show the comparison between the exact expression of the coefficient A 
and the expansions in the two regimes: threshold expansion on the left hand side and 
small-mass expansion on the right hand side. For the small-mass expansion we consider 
the limit m? / s — > for <f> = —{t — m 2 )/s = const and neglect terms suppressed by powers 
of m 2 /s larger than a given order. This prescription preserves the symmetry of the result 
under exchange of forward and backward directions, eft —> 1 — cf>, for a given order of the 
expansion. 

With the arXiv submission of this paper, we provide algebraic expressions and numer- 
ical implementations of the threshold expansion up to order (3 5 included and of the small 
mass expansion up to order (m 2 /s) 2 included. 



6. Conclusions and Outlook 



In this paper, we presented the analytic calculation of the two-loop leading color corrections 
to the heavy-quark production matrix element in the gluon fusion channel. The diagrams 
required to calculate this coefficient are all planar. The result presented here retains the 
exact heavy-quark mass dependence; no assumptions on the hierarchy between the mass 
scales involved in the problem was made. The formula for the coefficient A in Eq. fl2.G|), 
which was obtained in this work, was validated against analytic results valid in the small- 
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mass region [23]. In addition, we found numerical agreement with the exact analytic 
expression for the IR poles which was obtained in [35] . 

Our result represents a gauge invariant sub-set of the full two-loop corrections to 
the partonic process gg — > ti. In order to complete the analytic calculation of the two- 
loop corrections, it is necessary to calculate all of the fermionic diagrams, and the non- 
planar gluonic diagrams. A large part of the non-leading color coefficients in Eq. ( |2.G| ) 
can be calculated within the same calculational framework employed here. However, it is 
known that some of the two-loop diagrams appearing in the gluon fusion channel cannot 
be expressed in terms the HPLs functional basis. For example, some of the diagrams 
with a closed heavy-quark loop involve a "sunrise" -type subtopology with three equal 
massive propagators and an external momentum which is not on the mass shell of the 
internal propagators. Such a three-propagator graph can be written only in terms of 
elliptic integrals [53]. The use of numerical [24] or seminumerical [54,55] methods could 
thus be unavoidable in the evaluation of these diagrams. 

In order to obtain NNLO predictions for the total ti production cross section and for 
differential distributions, it is necessary to combine the two-loop virtual corrections with the 
already available one-loop corrections to the ti+(l parton) process and with the tree- level 
matrix elements for the process ti+(2 partons) [14-16]. These diagrams with additional 
partons in the final state contribute to infrared-divergent configurations where up to two 
partons can become unresolved. Their implementation requires the application of a NNLO 
subtraction method. The methods presently available [17-19] do not provide yet the full 
counterterms for a 2 — > 2 hadronic process involving massive partons. Two subtraction 
methods for the NNLO calculation of the top-quark pair production cross section were 
outlined recently [20,21]. 

In the light of these advancements, the calculation of the full NNLO corrections to top 
quark pair production in hadronic collisions is gradually becoming feasible. 
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